    RowSums = function(x){
        if(is.matrix(x)) return(rowSums(x))
        return(sum(x))
    }

    
    PMaxEstRough = function(ratio, K, k.max, gop, rep = 500, seed){
        
        d1 = 2^(2*K + 1)
        
        set.seed(seed)
        a = matrix(rbinom((K+1)*rep,1,0.5),rep,); colnames(a) = c(0:K)
        l = matrix(rbinom((K+1)*rep,1,0.5),rep,); colnames(l) = c(0:K)
        
        if(length(dim(ratio)) == 3){
            k = GetKGen2(a,l,ratio,k.max,K)
        } else k = GetK(a,l,ratio,k.max,K)
        k[k>1] = 1
        
        if(sum(k==1) > 0 & length(dim(ratio)) == 3){
            k1.index = which(apply(k, 1, function(x) sum(x==1)>0))
            for(i in k1.index){
                k[i,k[i,]==1] = sample(k[i,k[i,]!=1], sum(k[i,]==1), replace=TRUE)
            }
        }
        
        if(sum(k==1) > 0 & length(dim(ratio)) == 2){
            k[k==1] = sample(k[k!=1], sum(k==1), replace = TRUE)
        }
        
        f = function(x, k, gop){
            # prod(k) * (x^d1) / prod(1-x*k) - gop
            sum(log(k)) * d1 / rep  + d1 * log(x) - 
                sum(log(1-x*k)) * (d1-1) / rep -
                log(1-x) - log(gop)
        }
        ff = function(x){
            # prod(k) * (x^d1) / prod(1-x*k) - gop
            RowSums(log(k)) * d1 / rep  + d1 * log(x) - RowSums(log(1-x*k)) * 
                (d1-1) / rep - log(1-x) - log(gop)
        }
        tol = 10
        if(length(dim(ratio)) == 3){
            p.max = rep(1-0.1^tol, dim(ratio)[1])
            temp  = ff(1-0.1^tol)
            non.trivial = which(temp > 0 & is.finite(temp))
            for(i in non.trivial){
                p.max[i] = uniroot(f, k = k[i,], gop = gop[i], lower=0, upper=1, 
                                   tol=0.01 / rep)$root    
            }   
            p.max[temp==Inf] = 0
        }else{
            p.max = 1
            if(ff(1-0.1^tol) > 0) 
                p.max = uniroot(f, k = k, gop = gop[i], lower=0, upper=1, 
                                          tol=0.01 / rep)$root    
        }
        
        return(p.max)
        
    }
    
  
    Judge = function(ratio, K, k.max, gop, rep = 500){
        
        d1 = 2^(2*K + 1)
        
        # set.seed(0)
        a = matrix(rbinom((K+1)*rep,1,0.5),rep,); colnames(a) = c(1:K,0)
        l = matrix(rbinom(K*rep,1,0.5),rep,);     colnames(l) = c(1:K)
        
        k = GetK(a,l,ratio,k.max, K)
        k[k>1] = 1
        
        f = function(x){
            # prod(k) * (x^d1) / prod(1-x*k) - gop
            sum(log(k)) * d1 / rep  + d1 * log(x) - sum(log(1-x*k)) * d1 / rep -
                log(gop)
        }
        
        f(1-min(1e-10,0.001/d1))
        
    }
    
    
    
    